RegPT: Direct and fast calculation of regularized cosmological power spectrum at 

two-loop order 
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We present a specific prescription for the calculation of cosmological power spectra, exploited here 
at two-loop order in perturbation theory (PT), based on the multi-point propagator expansion. In 
this approach power spectra are constructed from the regularized expressions of the propagators that 
reproduce both the resummed behavior in the high-fc limit and the standard PT results at low-fc. 
With the help of iV-body simulations, we show that such a construction gives robust and accurate 
predictions for both the density power spectrum and the correlation function at percent-level in 
the weakly non-linear regime. We then present an algorithm that allows accelerated evaluations of 
all the required diagrams by reducing the computational tasks to one-dimensional integrals. This 
is achieved by means of pre-computed kernel sets defined for appropriately chosen fiducial models. 
The computational time for two-loop results is then reduced from a few minutes, with the direct 
method, to a few seconds with the fast one. The robustness and applicability of this method are 
tested against the power spectrum cosmic emulator from which a wide variety of cosmological 
models can be explored. The fortran program with which direct and fast calculations of power 
spectra can be done, RegPT, is publicly released as part of this paper. 

PACS numbers: 98.80.-k, 98.65. Dx 



I. INTRODUCTION 



Since recombination, the large-scale structure of the 
Universe has evolved dominantly under the influence of 
both the cosmic expansion and the force of gravity act- 
ing on a pressure-less fluid. The statistical nature of its 
spatial clustering is therefore expected to bring valuable 
cosmological information about the dynamics of the cos- 
mic expansion and structure formation. Of particular 
importance is the measurement of baryon acoustic oscil- 
lations (BAOs) imprinted on the power spectrum or two- 
point correlation function (e.g., [IHEj]) from which one 
can precisely determine the cosmological distance to the 
high-redshift universe, and henceforth clarify the nature 
of late-time cosmic acceleration (e.g., [6l-[iq| ) . Precious 
information regarding the growth of structure are and 
will also be obtained from redshift-space distortions (e.g., 
[HHIijl) and weak lensing measurements (see [T(|[l7| and 
review papers [H, at scales ranging to the linear or 
quasi-linear to the non-linear regimes. This could be cap- 
tured with unprecedented details with the ongoing and 
future surveys, thanks to their redshift depth and large 
angular area, such as the Sloan Digital Sky Survey III 1 , 



the WiggleZ survey 2 , the Subaru Measurement of Imag- 
ing and Redsfhits 3 , the Dark Energy Survey 4 , the Big- 
BOSS project 5 , the Physics of the Accelerating Universe 
collaboration 6 and the ES A/Euclid survey 7 . 

With the advent of such wealth of observations, there 
is therefore a growing interest in the development of the- 
oretical tools to accurately compute the statistical quan- 
tities of the large-scale structure. At decreasing redshift 
and scale, the evolution of the large-scale structure how- 
ever deviates significantly from the linear theory predic- 
tion and non-linear gravitational clustering effects have 
to be taken into account. While A-body simulations can 
be relied upon in specific cases, because of the range of 
scales to be covered and the variety of models to ex- 
plore, they should be complemented by analytical inves- 
tigations that aim at computing the statistical properties 
of the large-scale structure from first principles, hence- 
forth extending the validity range of linear calculations. 
It is to be noted that even at the scale of BAOs, linear 
calculations and one-loop standard PT corrections per- 



1 www.sdss3.org 



2 wigglez.swin.edu.au 
sumire.ipmu.jp/en/ 

4 www . darkenergysurvey . org 

5 bigboss.lbl.gov/index.html 

6 www . pausurvey . org 

7 www.euclid-ec.org 
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form poorly (see e.g., [20I4221 ] ) asking for more advanced 
PT calculations. The improvement of perturbation the- 
ory is thus a critical issue for the scientific exploitation of 
the coming surveys. Various resummation schemes have 
been proposed in Refs. [20l. |22h33| | that aim at improving 
upon standard schemes. The aim of this paper is not to 
compare them but to propose, and test, a specific scheme 
that can be used routinely in practice. 

In this paper, we are particularly interested in one of 
the resummation treatments, advocated in Ref. [34| . In 
this approach, the standard PT expansion is re-organized 
by introducing the multi-point propagators. These are 
the ensemble average of the infinitesimal variation of the 
cosmic fields with respect to the initial conditions. A key 
property shown in the previous reference is that all the 
statistical quantities such as power spectra and bispec- 
tra can be re-constructed by an expansion series written 
solely in terms of the multi-point propagators. This is 
referred to as the multi-point propagator expansion or 
T-expansion. The advantage of this approach is that the 
non-pcrturbative properties, which can be obtained in 
standard PT by summing up infinite series of PT expan- 
sions, are whole encapsulated in the multi-point prop- 
agators, including the effect of vertex renormalization. 
Furthermore, the T-expansion has been found to be valid 
not only for Gaussian initial conditions, but also for non- 
Gaussian ones [35j. The construction of accurate calcu- 
lation scheme for power spectra and bispectra can then 
be split in pieces that can be tested separately. 

The second key property that leads us to consider such 
objects is that their global shape, e.g. their whole k- 
dependence, can be computed in a perturbation theory 
context and compared to iV-body results thanks to the 
high-fc exponential damping tail they all exhibit [34l.l36l|. 
All these properties make the multi-point propagators 
the most important building blocks in the T-expansion 
and the focus of our modeling efforts. In the following we 
will in particular make full use of the novel regularization 
scheme proposed in Ref. [37j that allows to consistently 
interpolate between standard PT results at low-fc and the 
expected resummed behavior at high-fc. This scheme has 
been explicitly tested for the two-point propagators up 
to two-loop order in Ref. [38| and for (specific shapes of) 
the three-point propagators in Ref. [371 ] . 

The first objective of this paper is to present an explicit 
calculation of the non-linear power spectrum and corre- 
lation function of the cosmic density field based on this 
regularized treatment. Of particular interest is the extent 
to which the proposed scheme for T-expansion works be- 
yond standard PT when corrections at next-to-next-to- 
leading, i.e. two-loop, order are included. Results will be 
checked with V-body simulations. We will see that the 
T-expansion with the regularized treatment of propaga- 
tors, which we hereafter call RegPT, has good convergence 
properties and agree remarkably well with simulations 
entirely covering the scales of BAOs at any redshift. 

The second objective of this paper is to design and 
exploit a method to accelerate the power spectrum com- 



putations. Power spectra calculations in the context of 
RegPT calculations are rather involved requiring multi- 
dimensional integrations that have to be done with time- 
consuming Monte Carlo calculations. Typically, com- 
puting the power spectrum at percent level from our 
scheme takes several minutes. While this is acceptable 
when a handful of models have to be computed, this 
is an obstacle when a large domain of parameter space 
has to be systematically explored. Making use of the 
T-expansion functional form, we found though that it is 
possible to exploit a novel technique for accelerated cal- 
culation, in which only one-dimensional integrals need to 
be evaluated while ensuring the same precision as rigor- 
ous RegPT calculations. The bottom line of this approach 
is to see the resulting nonlinear power spectrum as a 
functional of the linear power spectrum and then Tay- 
lor expand this form with respect to the linear spectrum 
shape. We found that for well chosen fiducial models, it 
is sufficient to Taylor expand to first order only. We are 
then led to prepare in advance a set of kernel functions 
encoding the RegPT results for well chosen fiducial mod- 
els, whose normalizations are left floating, from which 
the RegPT predictions for the target model can be calcu- 
lated. We publicly release the fortran code, RegPT, as a 
part of this paper 8 . 



The organization of this paper is as follows. We be- 
gin by recalling the basic equations for cosmic fluid and 
perturbation theory in Sec. [TXJ We introduce the multi- 
point propagator and give the power spectrum expression 
based on the T-expansion. With the regularized treat- 
ment of multi-point propagators, in Sec. IIII1 we examine 
the the power spectrum calculations including the correc- 
tions up to the two-loop order, and investigate their UV 
and IR sensitivity in evaluating the PT kernels. Then, in 
Sec. HVl a detailed comparison between PT calculation 
and TV-body simulation is presented, and the accuracy 
and range of validity of PT calculation is checked. Based 
on this, Sec. |V] describes in detail the method to acceler- 
ate the power spectrum calculations. Robustness and ap- 
plicability of the accelerated RegPT calculations to a wide 
range of cosmological models are tested against power 
spectrum cosmic emulator code in Sec. I VII Finally, in 
Sec. lVIIi we conclude and explore practical extensions of 
this work. The description of the publicly released code, 
RegPT, is presented in Appendix [Bl 



The code is available at 

www-utap.phys . s .u-tokyo . ac . jp/~ ataruya/regpt_code .html 
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II. EQUATIONS OF MOTION AND THE Ec 
T-EXPANSION [39 



([2]) can be recast as the integral equation (e.g., [2C 



) 



A. Equations of motion 

In what follows, we consider the evolution of cold 
dark matter (CDM) plus baryon systems neglecting the 
tiny fraction of (massive) neutrinos. Owing to the 
single-stream approximation of the collisionless Boltz- 
mann equation, which is thought to be quite accurate an 
approximation on large scales, the evolution of the CDM 
plus baryon system can be treated as an irrotational and 
pressure-less fluid system whose governing equations are 
continuity and Euler e quat ions in addition to the Pois- 
son equation (see Ref. [33] for review). In the Fourier 
representation, these equations are further reduced to a 
more compact form. Let us introduce the two-component 
multiplet (e.g., [H): 



f a (M) = (<5(fc;t), 



9(k;t) 

'~fW 



(1) 



where the subscript a = 1, 2 selects the density and the 
velocity components of CDM plus baryons, with S and 
6(x) = V-v(x) I (a H), where a and H are the scale factor 
of the Universe and the Hubble parameter, respectively. 
The function f(t) is given by f(t) = d\nD(t)/dlna, and 
the quantity D(t) being the linear growth factor. Then, 
in terms of the new time variable r\ = \nD(t), the evolu- 
tion equation for the vector quantity ^ a (k;t) becomes 



d_ 

J ab o 

or] 



d fci d k 2 ki k 2 )jabc(ki,k 2 ) 



(2tt 

x # 6 (fci;j7) ^ c {k 2 ;v) 



(2) 



where we used the summation convention, that is the rep- 
etition of the same subscripts indicates the sum over the 
whole multiplet components. In the above, the quantity 
8d is the Dirac delta function, and the time-dependent 
matrix Q a i,(r}) is given by 







1 



ttab(v) 



(3) 



with the quantity Q m (r]) being the density parameter of 
CDM plus baryons at a given time. The vertex function 
jabc becomes 



Jabc(kl,k 2 ) = < 



+ 5 ( fl ' 6 > c ) = (i» 2 >i) 



(4) 



2|fci| 2 |fe 2 | 2 





(o,6,c) = (2,2,2) 



otherwise. 



*a(fc; rj) = g a b(ri, r/o) <f>b(k) 
d 3 kx d 3 k 2 



5 D (k - fci - fc 2 ) Jbcd(ki,k 2 ) 



(2tt)3 

x # c (fci;r/)* d (fc 2 ;?/)- 



(5) 



The quantity <f> a {k) = Sl^affc, 770) denotes the initial con- 
dition, and the function g a b denotes the linear propagator 
satisfying the following equation, 



Sal 



d_ 

dr) 



Qabiv) 



gbc(v,v') = 0, 



(6) 



with the boundary condition gabiViV) — ^ab- The statis- 
tical properties of the field ^ a are encoded in the initial 
field (f> a , for which we assume Gaussian statistics. The 
power spectrum of <j) a is defined as 

(Mk)Mk')) = (2tt) 3 S D (k + k') P ab ,o{k). (7) 

In what follows, most of the calculations will be made 
assuming the contribution of decaying modes of lin- 
ear perturbation can be neglected. This implies that 
the field 4> a (k) is factorized as <j> a {k) = So(k)u a with 
u a = (1,1), and thus the initial power spectrum is writ- 
ten as P a b,o( k ) = Po(k)u a u b . 

Using the formal expression (|S|), a perturbative solu- 
tion is obtained by expanding the fields in terms of the 
initial fields: 



(8) 



The expression of the solution at each order is written as 

*i n Hk;v)= I ^^ Mfc-*i *n) 



., (2 7 r) 3 (™- 1 ) 
xji n \k 1 ,k 2 ,--- ,fe„;tj)5 (fci)---«o(fc»). (9) 

The kernel J 7 ^ is generally a complicated time- 
dependent function, but can be constructed in terms of 
the quantities jabc and g a b- Examples of the solutions are 
shown diagrammatically in Fig. [TJ Because we are inter- 
ested in the late-time evolution of large-scale structure 
only, we can take the limit rjo —00. As a consequence, 
the fastest growing term is the only surviving one and 
the kernel is simplified into 

J<")(fc 1 ,--- ) fc„) = e»^(») ym (fc 1) -.-,fc n ), (10) 

(n) 

where the function i^sym is the symmetrized standard 

PT kernel, sometimes written as i^i'sym = (F n , G n ), 
whose explicit expressions are obtained from recursion 
relations as recalled in I39I1. 



4 



— o 



where we introduced the shorthand notation, 



Sab 01) 



YabcOo, k2) 




FIG. 1: Diagrammatic representation of the standard PT ex- 
pansion. 



B. T-expansion and regularized PT treatment 

In this paper, we are more specifically interested in the 
power spectra P a b(k; rj), defined as 



9 a (k;rj)* b (k';t})) = (2irf 5 D {k + k') P ab (\k\; v ) . (11) 



Substituting a set of perturbative solutions ([S]) into the 
above definition, it is straightforward to obtain the suc- 
cessive perturbative expressions for the power spectra. 
This is the standard PT treatment where the initial fields 
values are seen as the perturbative variables. The stan- 
dard PT calculations have, however, been shown to pro- 
duce ill-behaved higher-order corrections that lack good 
convergence properties. 

As an alternative to the standard PT framework, it 
has been recently advocated by many authors that the 
PT expansion can be re-organized by introducing non- 
perturbative quantities to improve the resulting conver- 
gence of the expansion. The T-expansion is one such 
non-perturbative framework, and the so-called multi- 
point propagators constitute the building blocks of this 
T-expansion. Denoting the (p + l)-point propagator by 
r( p ', we define 



1 



(fei) • • • S^ikj, 



Su(k - fei... p ) 



(2tt) 3 (p- 1 ) 



c p (ki 7 • • ■ , k p ; rj). 



(12) 



With these objects, the power spectra is shown to be 
expressed as [34j, 

P ab (| fc ];r ? ) = gt! / ^ffi Mfc-<h-,) 

X rW( Ql , • • • , q t ; ??)lf • • • , q t - V ) P (<?l) • • ' Po(Qt), 

(13) 



rW( Ql ,---, 9t ;r ? ) = r 



(*) 



,q t ;r))u ci ■■■u c 



(14) 

The diagrammatic representation for multi-point propa- 
gator and the power spectrum is respectively shown in 
Figs. [U and 

The construction of the T-expansion is rather trans- 
parent, and like Eq. (IT3l) . one easily finds the expressions 
for the higher-order statistical quantities such as bispec- 
trum. Another important point is that one can exploit 
the asymptotic properties of the propagators r( p ) beyond 
perturbation theory expansions. To be precise, in the 
high-fc limit, higher-order contributions can be systemat- 
ically computed at all orders, and as a result of summing 
up all the contributions, the multi-point propagators are 
shown to be exponentially suppressed [34J, |36[ , 



ri p) (fci, 



• , k p ; rf) 
exp 



,(p) 



,fc p ;r?) (15) 



with k = \k\ ^ \-k p \. This is the generalization of the 

result for the two-point propagator in Ref. [23j . Here, the 
quantity r^ reo is the lowest-order non- vanishing propa- 
gator obtained from the standard PT calculation, and o"d 
is the one-dimensional root-mcan-square of the displace- 
ment field defined by, 



a* 1 
ffd= 3 



cPq P (q) 



(2tt) 3 q 2 



(16) 



The form of Eq. (IT51) does not however provide a good 
description of the propagators at all scale. At low-fc the 
propagators are expected to approach their standard PT 
expressions that can be written formally, 

rW(fe ir - ) fe p ;r,) = rW ee (fe 1) --. ) fe p ;r,) 

oo 

+ E r ilioop(fci>---,^). (17) 

n=l 

For the dominant growing-mode contribution we are in- 
terested in, each correction term is expressed in terms of 
the standard PT kernels as, 



" (p) (k, ■ ■ ■ 
a, tree ^Ij 



,k p -,r,)=e^F^){k 1 ,c...,k p ), (18) 



for the tree-level contribution, and 
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FIG. 2: Example of the multi-point propagator, r„ . A large filled circle symbolically represents all possible contributions 
that enter into the fully non-linear propagator. A part of those contributions can be seen graphically using PT expansion (see 
Figs. [4] and [5] for three-point propagator ri 2 '). 




FIG. 3: Diagrammatic representation of the power spectrum by means of T-expansion. Here, the result up to the two-loop 
order is shown. In each contribution of the diagrams, the multi-point propagators are "glued" together at the crossed circles 
where the initial power spectra Po(fc) are inserted. 



rS- lo op(fcl: ~',k p ;v) = e (2n+P) " 4 P) / ^I^yI^ F^ytFHPl, 'Pi,' - ,P n , ~P n M, k p )P ( Pl ) ■ • ■ P (p n ) 

^ e ( 2n+P ), r (p)_ ioop(A;ij ... jM (1Q) 



for the n-loop order contributions, where the coefficient 
cfi 1 is given by = (2n+ P )Cp (2n — 1)!! where is the 
binomial coefficient. The graphical representation of the 
standard PT expansion is shown in Fig. [J] The important 
remark in Eq. (|19p is that each perturbative correction 
possesses the following asymptotic form, 

r ( P ) fc^oo 1 / k 2 ale 2r > \ n (p) 

a,n-loop | I 2 / a,n-treei V^"^ 

which consistently recovers the expression (|15[) when we 
sum up all the loop contributions. This indicates the exis- 

I 



which consistently reproduces one-loop PT results at low- 



tence of a matching scheme which smoothly interpolates 
between the low-fc and high-fc results for any multi-point 
propagator. Such a scheme has been proposed in Ref [37j 
where a novel regularized scheme, in which the low- and 
high-fc behaviors are jointly reproduced, is derived. The 
construction of the regularized propagator is totally un- 
ambiguous. They can incorporate an arbitrary number 
of loop corrections. 

Restricting the results to the growing mode contri- 
butions, the regularized propagators are expressed in a 
transparent way in terms of the standard PT results, and 
one gets 



k. An example of the regularized propagator valid at one- 



v n 



a,sym 



(fc 1 ,...,fc I ,)|l + ^^!|+e 2 Tjl. loop (fc 1> ... ,fc p ) expj 



* ^ ( (21) 
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loop order is diagrammatically shown in Fig.O This con- propagator including the corrections up to the two-loop 
struction is easily generalized to include the higher-order order becomes, 
PT corrections at low-fc. For instance, the regularized 



r i p L( fc ir-- ,k p ;rj) = e 



p. V 



F 



(v) 



k 2 a 2 e 2 ^ 



k 2 a 2 e 2r > 



' e2v rl-i oop (fci, • • • , k p ) \ 1 



k 2 a 2 e 2 ^ 



e i 2-loopl K l: 



exp 



(22) 



Note that the functions r„_ loop are the scale-dependent 
part of the propagator defined by Eq. p^|) . 



III. POWER SPECTRUM CALCULATION 
FROM REGULARIZED T-EXPANSION 

A. Power spectrum at two-loop order 

Since the proposed regularized propagators preserve 
the expected low-fc and high-fc behaviors, the convergence 



of the T-expansion adopting the regularization scheme 
would be much better than the standard PT expansion. 
In this paper, applying this regularized PT treatment, we 
will explicitly demonstrate the power spectrum calcula- 
tions at two-loop order, and comparing those predictions 
with iV-body simulations, the validity and precision of 
PT treatment are discussed. From Eq. (flU)) , the explicit 
expression for the power spectrum valid up to the two- 
loop order becomes, 



P afe (fc;r ? ) = rW eg (fc;r ? )rg eg (fc; ?7 )Po(A;)+2 J ^ 2 } cg (q, k - q; v)tfl s (q, k - q; V )P (q)P (\k - q\) 

+ 6 J T { a % g (p, q,k-p-q; r?)rg eg (p, q,k p - q; V )P Q (p)P Q (q)P (\k - p - q\) (23) 

with the regularized propagators given by 



r« (k;v) = e r > 



m c Jq,k-q; V ) 



k 2 a 2 e 2 ^ 



e 4 "r (1) 



a,2-loop 



(k) 



exp 



i k 2 a 2 e 2ri 
FW(q,k-q)ll + —£— 



e 2r Ti 2 ' (q,k-q) 



exp 



k 2 a\e 2r i 



r i 3 re K (P 5 Q,k-p- q\r)) = e 3v F$ (p, q,k - p - q) exp 



k 2 a 2 e 2 ^ 



k 2 a 2 e 2r > ' 



(24) 
(25) 

(26) 



Note that the higher-order contributions up to the two- 
and one-loop order of the propagators are respectively 
included in the expression of the regularized propagators 
ri^eg and ri 2 reg, consistently with the T-expansion at 
two- loop order. 

The power spectrum expression involves many inte- 
grals, but, most of them are reduced to two- or three- 
dimensional integrals if one uses the analytic expressions 
for the kernels of higher-loop corrections T in the regu- 



larized propagator. We use the expression in Ref. [38[ to 

—(2) 

evaluate r^^p, and adopt the fitting functions for the 

kernel of r^^p (see [38j]). We then apply the method 
of Gaussian quadrature to the numerical evaluation of 
the low-dimensional integrals. A bit cumbersome is the 
integral containing T^ 3 K While it can be reduced to a 
four-dimensional integral in principle, the expression of 
the resulting kernel would be very cumbersome and might 
not be suited for practical calculation. We thus adopt the 



(2) k 

Ta (pi,p 2 ) : — «- 



+ 




+ 




FIG. 4: Diagrammatic representation of the standard PT expansion for three-point propagator, ri 2 '. For fastest growing-mode 
contribution, the standard PT kernels, F^ym, form the basic pieces of PT expansion, depicted as incoming lines connected to 



a single outgoing line at the shaded circle. In the case of ri 2 ', the leading-order contribution is F^iym, and successively the 
kernels Fa%m and P<wm appear as higher-order contributions, for which pairs of the incoming lines are glued at the crossed 
circle, which indicates the initial power spectrum Pq, forming closed loops. 



(2) k 
Ta,reg (Pi, Pat : — «" 



X 



1 



-(kad)/2 

-O 



+ 




X 



i+O + tCO + t.OOO 



FIG. 5: Diagrammatic representation of the regularized three-point propagator, rii.reg- In the high-fc limit, the higher-loop 
contribution for three-point propagator behaves like Eq. (|20p . indicating that each loop diagram (r^_ loo ) is effectively split into 

tree diagram (Pifiym) and self-loop diagram ({— (ka^) 2 /2} n /n\), the latter of which is depicted as open loops. Systematically 
summing up all the higher-loop contributions, we recover Eq. (|15[l , which is graphically shown as the infinite sum of open-loop 
diagrams in the brace. To reproduce the standard PT result at low-fc, the one-loop diagram is inserted in the bracket, and the 
tree diagram is multiplied by the counter term {1 + (fcffd) 2 /2}. 



Monte Carlo technique of quasi-random sampling using 
the CUBA library [4fJ , and evaluate the five-dimensional 
integral directly 9 . 

Fig. |5] illustrates an example how each correction term 
in the regularized T-expansion contributes to the total 
power spectrum. The plotted result is the density power 
spectrum, P 11; and the contribution of the term involv- 
ing each multi-point propagator is separately shown. The 
three corrections contribute to the power spectrum at dif- 

(2) 

ferent scales, and the higher-order terms involving T rC g 

(3) 

and Treg are well-localized, each producing one bump. 
This is a clear manifestation of the result of the regu- 
larized PT treatment, and it resembles what the RPT 
calculations by Ref. [2J| give. 

In the next section, the results of the regularized T- 
expansion will be compared with TV-body simulations. 
But, before doing that, we will give several remarks and 
comments on the computation of the power spectrum in 
the subsequent subsection. 



Since the final result of the integration is expressed as a function 
of only the wavenumber k, the integrand possesses an azimuthal 
symmetry with respect to the vector, k, indicating that the in- 
tegral is reduced to a five-dimensional integral. 



B. Effect of running UV cutoff for era 

Since the shape of the power spectrum given by 
Eq. (|23j) significantly depends on the exponential damp- 
ing in the regularized propagators, we first comment on 
the effect of this function. As it has been shown, the ex- 
ponential function arises from the summation of infinite 
series of perturbations at all order in the high-fc limit. 
Recently, Ref. [36[ advocated that this exponential func- 
tion can be interpreted as the result of resummation at 
hard part (high-fc), and the displacement dispersion 
in the exponent must be evaluated in a consistent way 
that the domain of the integral is restricted to a soft part 
(low-fc). This implies that depending on the scale of our 
interest, the boundary of the soft and hard domains can 
be changed, and the resulting quantity Od should be re- 
garded as a scale-dependent function. 

In Fig. [JJ we examine the impact of the scale- 
dependent <7<j on the power spectrum at z = 1. Plotted 
results are the contributions of the power spectrum cor- 
rections (upper) and the total power spectrum divided by 
the smooth reference linear spectrum (bottom) . Here, we 
evaluate by introducing the running UV cutoff k\(k): 




Various curves in Fig. [7J represent the results with dif- 
ferent prescription for the running UV cutoff. The cor- 
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FIG. 6: Contribution of multi-point propagators to the power spectrum, P(k) — Pn(k) at z = 1. Magenta, green, and blue 
curves represent the power spectrum contributions from the first, second, and third terms at the right-hand-side of Eq. (|23jl . 

respectively, each of which just corresponds to the diagram in Fig. [3l involving r[lg, Treg, and r^g. Summing up these 
contributions, total power spectrum is shown in black solid line. For reference, linear power spectrum is also plotted as dotted 
line. 



(3) 

rcction involving the four-point propagator Trcg is most 
sensitively affected by the running cutoff, and the re- 
sulting power spectrum significantly varies at scales k > 
0.2 ZiMpc -1 . This is because the exponential damping 
manifests itself at the scale k ~ l/(e v ad) where the con- 

(3) 

tribution from If-cg, which contains no relevant terms 
counteracted with the exponential damping, becomes sig- 
nificant among the three corrections. 

In the bottom panel of Fig. [7J we also plot the result 
of TV-body simulations (see Sec. II V A"|) . The comparison 
with simulation suggests that the PT calculation with 
running cutoff k\ ~ k/2 — fc/5 is favored, although there 
is no clear physical reason why this is so. Strictly speak- 
ing, the running IR cutoff might also be introduced in 
evaluating all the integrals in the power spectrum ex- 
pression, so as to consistently discriminate between the 
contributions coming from soft and hard parts. More- 
over, the running cutoff k\ may also depend on the red- 
shift. These complications mostly come from the ambi- 
guity of the boundary between soft and hard domains in 
our regularization scheme. For practical purpose to the 
cosmological application, we postpone these issues to fu- 
ture investigation, and take a rather phenomenological 
approach. Hereafter, the running cutoff is only intro- 
duced in evaluating a^, and we evaluate it according to 
Eq. (JUJ, setting the cutoff scale to fc A = k/2. With 



this treatment, we will see later that the PT prediction 
becomes improved compared to the standard PT calcu- 
lation, and it reproduces the iV-body results quite well 
at any redshift. 



C. Sensitivity to IR and UV cutoff 

In computing the power spectrum, except for ad, the 
domain of each integral in Eq. (I2U1) are usually taken 
broad enough so as to ensure the convergence of the re- 
sults. In comparison with iV-body simulations, however, 
a care must be taken because the available Fourier modes 
in simulations are restricted depending on the simulation 
box size and/or mesh size of Fourier transform, which af- 
fect both the efficiency of mode transfer and strength of 
mode coupling. The evolved result of the power spectrum 
would thus be changed, and it should be carefully com- 
pared with PT calculation, taking the finite resolution 
into consideration. Here, focusing on the BAO scales, we 
briefly discuss the sensitivity of the PT calculation to the 
IR and UV cutoff in the integrals. 

Fig. [8] shows the variation of the power spectra with re- 
spect to the UV (left) and IR (right) cutoff. In general, 
the kernel of integrals becomes broader for higher-loop 
corrections, and thus the two-point propagator F^g con- 
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FIG. 7: Sensitivity of the power spectrum prediction at z = 1 to the UV cutoff in the estimation of era. Top panel shows each 
contribution of the power spectrum corrections involving I^ 1 ' (magenta), r' 2 ' (green), and I^ 3 ' (blue), respectively. Bottom 
panel shows the total sum of power spectrum divided by the smooth reference power spectrum, P no -wiggie(ft), which is calculated 
from the no- wiggle formula of the linear transfer function in Ref. In each case, top lines represent the results obtained 

by setting era = 0, while undermost lines show the cases adopting the value of era without UV cutoff. The middle six lines 
represent the cases adopting the running UV cutoff in estimating a^, with cutoff fcA(fc) = k, k/2, k/3, k/5, fc/10, and fe/20 
(from bottom to top). As a reference, linear theory prediction is also plotted in both panels (dotted). 



taining the two-loop contribution is sensibly affected by 
the UV cutoff. Note that the signs of one- and two-loop 
corrections in Treg are opposite at BAO scales. Hence, as 
decreasing the cutoff wavenumber £; max , the cancellation 
of each contribution is relatively relaxed, and the power 
spectrum amplitude gets increased. On the other hand, 
due to the lack of the long-wave modes, the IR cutoff not 
only decreases each contribution of the loop integrals, 
but also reduces Od, leading to a slight suppression of 
the exponential damping. The net effect of the IR cutoff, 
especially at small scales k > 0.2 h Mpc -1 , is that the lat- 
ter overcomes the former, and the total power spectrum 
is slightly enhanced. 

These results imply that the effect of UV and IR cut- 
off not only affects the power spectrum shape at small 
scales, but also causes a slight offset in power spectrum 
amplitude at moderately large scales, k ~ 0.1/iMpc -1 . 
The size of these effects is basically small, but would not 

I 



be negligible in a percent-level comparison. Based on 
this remark, in what follows, we adopt the cutoff scales 
(kmin, fcmax) = (5 x 10~ 4 , 10) h Mpc -1 as default param- 
eters to compute the power spectra. With this setup, 
RegPT calculation gives a mostly convergent result, which 
can be compared with high-resolution iV-body simula- 
tions with a large box size. 

D. Comparison with MPTbreeze 



In Ref. [42J, mptbreeze, an alternative scheme has 
been proposed for the construction of power spectra that 
is based on the same multi-point propagator expansion. 
This proposition is, however, based on simplified assump- 
tions regarding the behavior of the multi-point propaga- 
tors. More specifically, in mptbreeze, the propagators 
are assumed to take the form, 



ltL g (<Z l5 • ■ • , q p -x, k - q 1 ... (p . 1) ;v) = e pri F^% m ( qi , q p _ 1: k - q,...^) exp{f a (k)e 2 ^} , (28) 

I 

where fi(k) and f2(k) are the one- loop corrections to the corresponds to the late-time original expression of the ex- 
density and velocity propagators, respectively. This form 
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FIG. 8: Sensitivity of the power spectrum prediction to the UV (left) and IR (right) cutoff. In each panel, the ratio of power 
spectrum, P(k)/P no - vr i m i e (k), is plotted as function of k at z = 1 (left) and 0.35 (right). In evaluating the integrals of the 
power spectrum corrections, the maximum wavenumber for the range of the integral is set to 2-k (green), 7r (cyan), and 7r/2 
(blue) ftMpc -1 in the left panel, while in the right panel, we change the minimum wavenumber fc m i n to k m - ln — 27r/Lbox with 
ibox = 2,048 (green), 1,024 (cyan), and 512 h' 1 Mpc (blue). The magenta lines indicate the results adopting the default set 
of parameters (fc m in, & max ) = (5 x 10~ 4 , 10) h Mpc -1 . 



ponentiation scheme initially put forward in Ref. [23j . It 
is shown in Refs. [37, 38] that at one-loop order, this pre- 
scription gives nearly identical result for the two-point 
propagator to the prescription proposed in Ref. [23| . 
The mptbreeze prescription, however, ignores the im- 
pact of two-loop PT corrections on the two-point propa- 
gators. From the results presented in [38], it implies that 
mptbreeze might be outperformed by RegPT at z > 1. 
On the other hand, the predictions of that scheme are 
made more robust because they are less sensitive to the 
UV part of the linear spectrum as discussed in that paper. 
Furthermore, the one-loop correction for the three- and 
four-point propagators are treated in an effective way. 
These simplified assumptions allow a more rapid calcu- 
lations of the set of diagrams. It takes just a few seconds 
to get the expected shape in this scheme. The compu- 
tational time is, however, rather comparable to the fast 
implementation of RegPT which we will present in Sec.fVl 



IV. COMPARISON WITH TV-BODY 
SIMULATIONS 

We are now in position to present quantitative com- 
parisons between RegPT calculations and TV-body simu- 
lations. After briefly describing the TV-body simulations 
in Sec. IIV Al we show the results of power spectrum and 
two-point correlation function in Sec. IIV FJI and ITVC] re- 
spectively. Precision and validity of the PT predictions 
are discussed in detail. 



A. TV-body simulations 

To compare RegPT calculations with TV-body simula- 
tions, we ran a new set of TV-body simulations, which will 
be presented in more detail with an extensive convergence 
study in Q. This set of simulations can be regarded as 
an updated version of the one presented in [22] with much 
larger total volume and a more careful setup to achieve 
a smaller statistical and systematic error. The data were 
created by a public TV-body code GADGET2 [H| with cubic 
boxes of side length 2, 048 /i _1 Mpc, and 1, 024 3 particles. 
The cosmological parameters adopted in these TV-body 
simulations are basically the same as in the previous one, 
and are determined by the five-year WMAP observations 
[4rj| (see Table IJ) . The initial conditions were generated 
by a parallelized version of the 2LPT code [47l|. d eveloped 
in Ref. [48] . After several tests given in Ref. [48[ , a lower 
initial redshift z; n ;t turns out to give a more reliable esti- 
mate for the power spectrum at BAO scales, and we thus 
adopt the initial redshift Zj n i t = 15. With this setup, we 
have created 60 independent realizations and the data 
were stored at redshifts z = 3, 2, 1, and 0.35. The total 
volume at each output redshift is 515/i -3 Gpc 3 , which is 
statistically sufficient for a detailed comparison with PT 
calculations. 

We measure both the matter power spectrum and 
the correlation function. For the power spectrum, we 
adopt the Cloud-in-Cells interpolation, and construct the 
Fourier transform of the density field assigned on the 
1, 024 3 grids. As for the estimation of the two-point cor- 
relation function, we adopt the grid-based calculation 
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FIG. 9: Comparison of power spectrum results between TV-body simulations and RegPT calculations. In each panel, the results 
at z = 3, 2, 1, and 0.35 are shown (from top to bottom). Left panel shows the ratio of power spectrum to the smooth 
linear spectrum, P(fc)/P no - w i gg i c (fc), where the reference spectrum P n o-wig g ie(fc) is calculated from the no- wiggle formula of 
the linear transfer function in Ref. [4l| . Solid lines are the RegPT results, while dotted lines represent the linear theory 
predictions. Right panel plots the difference between A-body and RegPT results normalized by the no-wiggle spectrum, i.e., 
[^N-body(fc) — -PRcgPT(fc)]/Pno-wiggie(fc)- In each panel, The vertical arrows respectively indicate the maximum wavenumber 
below which a percent-level agreement with JV-body simulation is achieved with Lagrangian resummation theory [25l . |43| and 
closure theory [22J, |29(, including the PT correction up to two- loop order. 



TABLE I: Cosmological parameters for A-body simulations (ACDM) 



Name 


ibox 


# of particles 




# of runs 








w 


h 


n s as 


wmap5 


2,048fr~ 1 Mpc 


1, 024 3 


15 


60 


0.279 


0.721 


0.165 


-1 


0.701 


0.96 0.815 9 



using the Fast Fourier Transformation [22j . Similarly 
to the power spectrum analysis, we first compute the 
square of the density field on each grid point in Fourier 
space. Then, applying the inverse Fourier transforma- 
tion, we take the average over separation vectors and 
realizations, and finally obtain the two-point correlation 
function. The implementation of this method, together 
with a convergence test, is presented in more detail in 
Ref .[13. 



B. Power spectrum 

Let us first present the power spectrum results. 
Left panel of Fig. |H] shows the ratio of the power 
spectra, P(k)/ P no - w [ gg i c (k), while the right panel 
plots the fractional difference between Af-body simu- 



lations and PT calculations, defined by [PN-body(^) — 
-PpT(fc)]/P n o-wiggie(fc) (where P n o-wig g io(fc) is calculated 
from the no- wigg le formula of the linear transfer func- 
tion in Ref. [41|). Overall, the agreement between 
RegPT and A^-body simulations is remarkable at low-fc, 
and a percent-level agreement is achieved up to a cer- 
tain wavenumber. For a decreasing redshift, the non- 
linearities develop and the applicable range of PT calcu- 
lations inevitably becomes narrower, however, compared 
to the standard PT predictions, the RegPT result can re- 
produce the A^-body trend over an even wider range. In- 
deed, the range of agreement with A^-body simulations 
is rather comparable to other improved PTs including 
higher-order corrections, such as closure theory [22t |29j. 
and better than some of those predictions. For refer- 
ence, we compute the power spectra from closure and La- 
grangian resummation theory (LRT) [2^, |43| at two- loop 



12 



0.0025 
0.002 
0.0015 

o.ooi 

0.0005 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



Linear 

RegPT (2-loop)~ 




10 3 



10 2 



60 70 80 90 100 110 120 130 140 150 
r [h"'Mpc] 




40 60 80 100 
r [h _1 Mpc] 



FIG. 10: Comparison of two-point correlation function between TV-body and RegPT results at z = 3, 2, 1, and 0.35 (from 
bottom to top). In each panel, magenta solid, and black dotted lines represent the prediction from RegPT and linear theory 
calculations, respectively. Left panel focuses on the behavior around baryon acoustic peak in linear scales, while right panel 
shows the overall behavior in a wide range of separation in logarithmic scales. Note that in right panel, the resulting correlation 
function is multiplied by the cube of the separation for illustrative purpose. 



order, and estimate the range of a percent-level agree- 
ment with iV-body simulations, the results of which are 
respectively depicted as green and blue vertical arrows in 
right panel of Fig. GO Note that at z = 3, the range of 
agreement for closure theory exceeds the plotted range, 
and is not shown here. 

Although the RegPT treatment gives a very good per- 
formance comparable to or even better than other im- 
proved PTs, a closer look at Fig. [9] reveals several sub- 
percent discrepancies. 

• One is the low-fc behavior at z — 0.35, which ex- 
hibits a small discrepancy with iV-body simulation. 
Our investigations indicate that it is probably due 
to a poor convergence of standard PT expansion, 
since the low-fc behavior of regularized propagators 
heavily relies on the standard PT treatment. To be 
specific, the convergence of Ireg is the main source 

of this discrepancy. Indeed, if T r l g is computed 
at one-loop order only, the power spectrum is en- 
hanced, and then N-body results at low-fc lie in be- 
tween the two predictions. The impact of the high 
order PT corrections to the two-point propagator 
are specifically studied in a separate publication, 
381. 



results at mid-fc regime (fc ~ 0.2-0.3/iMpc^ 1 ). 
It is unlikely to be due to a poor convergence of 
standard PT expansion. We rather think that the 
performances of the iV-body simulations might be 
responsible for this (small) discrepancy. We have 
tested several runs with different resolutions, and 
found that the low-resolution simulation with a 
small number of particles tends to underestimate 
the power at high-z. Possible reason for this comes 
from the precision of force calculation around the 
intervening scales, where the Tree and Particle- 
Mesh algorithms are switched, and is mainly at- 
tributed to the inaccuracy of the Tree algorithm. 
Though the intervening scale is usually set at a 
sufficiently small scale, with a low-resolution simu- 
lation, it may affect the large-scale dynamics with 
noticeable effects at higher redshifts. Systematic 
studies on the convergence and resolution of N- 
body simulations will be reported elsewhere 



• Another discrepancy can be found in the high-z 
results, which temporally overshoot the iV-body 
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C. Correlation function 

We next consider the two-point correlation function, 
which can be computed from the power spectrum as 



dk k 2 . sm{kr) 
P(k) ■ 



2ir 2 



ki 



(29) 



In Fig. HOI left panel focuses on the behaviors around 
the baryon acoustic peak, while right panel shows the 
global shape of the two-point correlation function plot- 
ted in logarithmic scales, for which £(r) has been multi- 
plied by the cube of the separation. The RegPT results 
agree with TV-body simulations almost perfectly over the 
plotted scales. As it is known, the impact of non-linear 
clustering on the baryon acoustic peak is significant: the 
peak position becomes slightly shifted to a smaller scale, 
and the structure of the peak tends to be smeared as the 

redshift decreases (e.g., (Il, [H IH [H3 ) ■ The Re S PT cal " 
culation can describe not only the behavior around the 
baryon acoustic peak but also the small-scale behavior 
of the correlation function. Note that similar results 
are also obtained from other improved PT treatments 
such as closure and LRT. Although the RegPT predic- 
tions eventually deviate from simulations at small scales 
- the result at z — 0.35 indeed manifests the discrepancy 
below r ~ 30/i -1 Mpc - the actual range of agreement 
between RegPT and TV-body results is even wider than 
what is naively expected from the power spectrum re- 
sults. In fact, it has been recently advocated by several 
authors that with several improved PT treatments, the 
one-loop calculation is sufficient to accurately describe 
the two-point correlation function (e.g., [H,|4J,[5l|). We 
have checked that the RegPT treatment at one-loop or- 
der can give a satisfactory result close to the two-loop 
result, and the prediction including the two-loop correc- 
tions only slightly improves the agreement with iV-body 
simulations at small scales. This is good news for practi- 
cal purposes in the sense that we do not necessarily have 
to evaluate the multi-dimensional integrals for the accu- 
rate prediction of two-point correlation function in the 
weakly non-linear regime. Nevertheless, in this work, we 
keep the two-loop contributions in the computed contri- 
butions. The computational costs of the two-loop order 
will be addressed in the following with the development 
of a method for accelerated PT calculation at two-loop 
order. 



REGPT-FAST: ACCELERATED POWER 
SPECTRUM CALCULATION 



In this section, we present a method that allows ac- 
celerated calculations of the required diagrams of the 
two-loop order RegPT prescription. In principle, the 
power spectra calculations in the context of RegPT re- 
quire multi-dimensional integrations that cannot be done 
before-hand as they fully depend on the linear power 
spectra. It is however possible to obtain the required 
quantities much more rapidly provided we know the an- 
swer for a close enough model. 

The key point in this approach is to utilize the fact that 
the nonlinear RegPT power spectrum is a well-defined 
functional form of the linear power spectrum. Each of 
the diagram that has to be computed is of quadratic, 
cubic, etc. order with respect to the linear power spec- 
trum with a kernel that, although complicated, can be 
explicitly given. It is then easy to Taylor expand each of 
these terms with respect to the linear power spectrum. 
In principle one then just needs to prepare, in advance, 
a set of the RegPT results for some fiducial cosmologi- 
cal models, and then take the difference between fiducial 
and target initial power spectra for which we want to cal- 
culate the non-linear power spectrum. These differences 
involve only one-dimensional integrals at the first order 
in the Taylor expansion. 

In the following, we present the detail of the implemen- 
tation of this approach illustrating it with the one-loop 
calculation case. 



Power spectrum reconstruction from fiducial 
model 



While our final goal is to present the fast PT calcula- 
tion at two-loop order, in order to get insights into the 
implementation of this calculation, we consider the power 
spectrum at one-loop order. The complete expressions 
needed for the fast PT calculation at two-loop order, 
together with the prescription how to implement it, is 
presented in Appendix |A"1 

Compared to the expressions given in Eq. (|23|) . the 
power spectrum at one-loop order of the T-expansion re- 
duces to 



P ab (k-,r 1 )=r^ g (k;r 1 )ri% s (k-,r ] )Po(k) + 2 J ASL T^ Ies (q,k - q;r,)4% s (q,k - q;r 1 )Po(q)Po(\k - q\) (30) 



with the regularized propagators r^g and T^g valid at one-loop order being: 



r£> (k;r,) = e r > 



k 2 aW n 



+ ^T ( -\_ x (k) 



x exp 



p 2 2r, 

W l (31) 



r$eg(g, k-q; V) = e 2v Fi% m {q, k - q) 
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Note that the quantity T i_ loop is defined in Eq. (|T9| . 
and explicitly given by 

/t3 
j^F(% m (q,-q,k)P (q). (33) 

Thus, in Eq. (|30|1 . there apparently appear two contribu- 
tions which involve multi-dimensional integrals; r a l . loop 

in the regularized propagator ri^-cg, and the second 
term at the right-hand side. Although these contribu- 
tions are known to be further reduced to one-, and two- 
dimensional integrals (e.g., Refs [H, [2i| H3]), respec- 
tively, for the sake of this presentation we keep the ex- 
pressions as in their original form. 

As it has been mentioned earlier, the key idea of accel- 
erated calculation is to prepare a set of RegPT results for 
fiducial cosmological models. Let us denote the initial 
power spectrum for fiducial cosmology by Po,ftd(k). And 
we denote the initial spectrum for the target cosmologi- 
cal model, for which we want to calculate the non-linear 
power spectrum, by -Po,tar g et(fc)- For the moment, we as- 
sume that the difference between those spectra is small 
enough. Then, we may write 

0, target (k)=P 0M (k) + SP (k). (34) 
Hereafter, we focus on the power spectrum of density 



field, P\\, and drop the subscript. Substituting the above 
expression into Eqs. (j3"0)) -(f3"2" j) . the non-linear power spec- 
trum for the target model is symbolically written as 

-Ptargct(fc; V) = fun-pert Vi a d, target I Po, fid (k)\ 

+ Pcorr[k i r],ad,teir g ct;SP (k)}. (35) 

Here, the first term at the right hand side is the un- 
perturbed part of the one-loop power spectrum, which 
is nothing but the expression (I30[) adopting the ini- 
tial power spectrum for fiducial model, i"b,fid(^), but 
with the cosmological dependence of the time variable, 
given by rj — ln£>(z), being calculated from the target 
model. Also, the dispersion of displacement field, Cd, 
should be replaced with the one for the target model, 
i.e., Od, target = [/ dqP QMrgct (<?)/(67r 2 )] 1 / 2 . In each term 
of Eq. (l30|) . the scale and time dependence can be sep- 
arately treated, and thus the un-perturbed power spec- 
trum, Pun-perturbed, is evaluated algebraically by sum- 
ming up each contribution, for which we use the precom- 
puted data set in evaluating the scale-dependent func- 
tion. 

In Eq. ([331) . the contribution P corr includes the non- 
linear corrections originating from the differences of ini- 
tial power spectra between fiducial and target cosmolog- 
ical models. To first order in SPq, we have 



-Pcorr [fc, ?7, cr dM r gc uSPo(k)} = 2T£>(k; rj) ST^H (k; 7?) Po,fid (k)+ rW(fc;77) 6P (k) 



d 3 q 



rW (q, k - q;rjj\ P , nd (|fe - q\)6P (q) 



(36) 



In the above expression, The quantity ad appearing in the 

propagators Ireg and Treg should be evaluated with the 
linear power spectrum for the target cosmological model. 

The perturbed propagator tflieg is expressed as 

5riH(k- v ) = e^ST {1 \k)e~ k2 <^ e2 "^ 2 (37) 

where the kernel of integral in ST ^ is the same one as 
in Eq. (I33p . but we may rewrite it with 

ST {1 \k)= I ^-L^(q,k)6P (q) (38) 



with the kernel given by 

4%,fc) = 3 J ^L F [% m {q,~q,k). (39) 

Since the kernel U± only includes the PT kernel whose 
cosmological dependence is extremely weak, we can sepa- 
rately prepare the numerical data set for in advance 



10 . Then, we can use it to compute ST for arbitrary 
SPq, where the remaining integral to be evaluated is re- 
duced to a one-dimensional integral. 

Furthermore, the integral in the last term of Eq. (|36[) 
is rewritten with 



_fq_ 

(2tt) 3 

_ „-fc ; 



TM(q,k-q;ri) P 0M (\k - q\) SP (q) 

2 XW(q,k)5P Q (q) (40) 



■:A>I 



dq q 2 



2tt 2 



Indeed, the kernel Lj 1 ' is analytically known, and the explicit 



expression is given in, e.g., Refs. 1231. |29|, 137 
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with the function being 



x P MVk 2 -2kqfi + q 2 ), (41) 

where the variable /i is the directional cosine defined by 
fi = (k ■ q)/(kq). In deriving the above expression, we 

(2) (2) 

used the symmetric property of r r0 g, i.e., r rcg (fci, ks) — 
r^eg(fe2) ki). Since the quantity X^(q,k) can be com- 
puted in advance, all the integrals involving the power 
spectrum <5P are shown to be effectively reduced to one- 
dimensional integrals. In other words the only remaining 
task is to evaluate one-dimensional integrals, which can 
be done very efficiently. 

The practical implementation of this method makes 
use of another important property of the kernel functions. 
They indeed have a very simple dependence on a global 
rescaling of the power spectrum, Po,fid — > c Po,fid- It 
is then possible, without extra numerical computation, 
to choose the fiducial model among a continuous set of 
models. The model we choose, that is the normalization 
factor c we take, is such that the difference SPo{k) is as 
small as possible in the wave-modes of interest. As we 
will see in Sec. I VII it makes the use of this method very 
efficient. 

Note that although the treatment depicted here does 
not give much impact on the computational cost of the 
one-loop calculation, we will explicitly show in Appendix 
|A1 that at two-loop order the PT corrections involving 
multi-dimensional integrals can be similarly reduced to 
one-dimensional integrals. In the following, we denote 
RegPT-f ast the implementation of this approach at two- 
loop order. 



B. Performances 

Let us now illustrate the efficiency of the 
RegPT-fast expansion. Based on the expressions 
given in Appendix El we calculate the power spectrum 
and correlation function at two-loop order. We adopt 
the best-fit parameters determined by the third-year 
WMAP result [52] as the fiducial cosmological model 
from which we try to reproduce the RegPT results for 
the five-year WMAP cosmological model. Cosmological 
parameters for the fiducial model is listed in Table [TTJ 
Compared to the target model in Table HI the mass 
density parameter shows a 20% difference, and with 7% 
enhancement in the power spectrum normalization (erg), 
this leads to a 20-30% difference in the initial power 
spectrum. 

Fig. Ql] plots the results of the RegPT-fast calculation 
(blue) compared to the target RegPT calculation (ma- 
genta). We plot, for a specific redshift z = 1, the ratio of 
the power spectrum to the smooth reference spectrum, 
P(fc)/P no -wiggie(fc), and correlation function multiplied 
by the cube of separation, r 3 f (r), in left and right panels, 



respectively. The RegPT-fast results perfectly coincide 
with RegPT direct calculation, even outside the range of 
agreement with iV-body simulations. 

Note that the perfect match between RegPT and 
RegPT-fast results is due to a large extent to the con- 
tributions of the higher-order PT in the correction, P cor r 
or £corr- This appears clearly in the plots of the linear 
theory correction, SPq = Po. target — Po,fid and its Fourier 
counterpart S^q (cyan long-dashed). As shown in cyan 
solid lines, the total contribution, i.e., the combination 
of the un-perturbed part plus linear theory correction, 
somehow resembles the result with direct RegPT calcula- 
tion, but exhibits a rather prominent oscillatory feature 
with slightly different phase in power spectrum, leading 
to a non-negligible discrepancy. Accordingly, in correla- 
tion function, the acoustic peak becomes enhanced, and 
the position of peak is shifted to a small separation. Note 
finally, that these results could only be achieved with the 
help of the rescaling properties of the kernel functions. In 
this particular case the fiducial model has been rescaled 
as Po,fid — >• 1.082 Po. fid- Rescaling is a key feature of the 
RegPT-fast method. It will be further discussed in the 
next section. 



VI. TESTING REGPT TREATMENT FOR 
VARYING COSMOLOGICAL MODELS 



The purpose of this section is twofold. Our first 
goal is to explore the validity and applicability of 
the RegPT-fast scheme. Having shown that the 
RegPT-fast approach can be used in one specific exam- 
ple, we now want to discuss the usefulness of this treat- 
ment from a more practical point of view. To be precise, 
we want to know how well the RegPT-fast treatment 
can reproduce rigorous RegPT calculation in a variety of 
cosmological models. 

Our second and natural goal is to test the 
RegPT scheme itself, whether from direct or fast calcu- 
lations, against N-body based predictions such that the 
cosmic emulator 11 . 

To do that, we have selected the 38 cosmological mod- 
els investigated in Ref. [53] for which we can use the 
publicly released code, cosmic emulator, that provides 
interpolated power spectra derived from iV-body simula- 
tions. Let us remind that the cosmological models con- 
sidered there are sampled from a wide parameter space 



11 http: //www. lanl . gov/projects/cosmology/CosmicEmu/ 
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FIG. 11: Example of the performances of the RegPT-fast approach compared to direct RegPT calculation. Left panel shows the 
power spectrum divided by the smooth reference spectrum, P(fc)/Pno-wig g ic(fc), while right panel plots the correlation function 
multiplied by the cube of separation, i.e., r 3 £(r). In both panels, the results at z = 1 are shown, together with TV-body 
simulations. The RegPT-fast results, computed with prepared data set for fiducial cosmological model, are plotted as magenta 
solid lines, which almost coincide with those obtained from the rigorous RegPT calculation (magenta solid). As shown in Eq. 1351 
the RegPT-fast results are divided into two contributions; un-perturbed part (Pun-pert or £ un -pert) adopting the wmap3 model 
as fiducial cosmology, and the correction part (P CO rr or £ CO rr) evaluated with the power spectrum difference (5Po. These are 
respectively plotted as green dashed and cyan dashed lines. For reference, we also show the linearly evolved result of power 
spectrum difference <5Po (cyan dashed) and the sum of the contributions P U n- pC rt + <5Po (cyan solid) in left panel, and their 
Fourier counterparts in right panel. 

TABLE II: Cosmological parameters for fiducial models used for the RegPT-fast calculation (see Sec. I VI A")) . 



Name 








w 


h 


n s 




wmap3 


0.234 


0.766 


0.175 


-1 


0.734 


0.961 


0.760 


M001 


0.430 7 


0.5692 


0.150 


-0.816 


0.597 7 


0.946s 


O.8I61 


M023 


O.I6O2 


0.839s 


O.I8I7 


-1.261 


0.8694 


0.901 6 


O.6664 



for fiat w-CDM cosmology, and lie within the range, 



A. Convergence of RegPT-fast 



0.120 < n m h 2 < 0.155, 
0.0215 < Q h h 2 < 0.0235, 
0.85 < n 3 < 1.05, 
-1.30 < w < -0.70, 
0.616 < cr 8 < 0.9. 



The concrete values of the cosmological parameters in 
each model are not shown here. Readers can find them 
in Table 1 of Ref . [53[ . Fig. [T2] shows the linear power 
spectra Po(fc) for the 38 cosmological models, which have 
been all produced with CMB Boltzmann code, camb [54| . 



Let us first examine the convergence of the power spec- 
trum calculations between RegPT and RegPT-fast treat- 
ments. We ran both the RegPT-fast and RegPT codes, 
and evaluate the fractional difference between these 
power spectra, defined by pRegPTfast(^)/pRegPT(fc) - 1- 
Collecting the results at z = 1 in each cosmological 
model, the convergence of the power spectrum calcula- 
tions for 38 models is summarized in Fig. Q21 In left 
panel we show the result when only one fiducial model, 
wmap3, is used. In that case RegPT-fast results tend to 
underestimate the results from rigorous RegPT calcula- 
tions at increasing k, and most of them eventually ex- 
ceed the 1% difference, indicated by the green dashed 
line. This is because the shape of the initial power spec- 
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trum in each target model is rather different from that 
in the fiducial model, and even adjusting the re-scaling 
parameter c cannot compensate a large power spectrum 
difference. To be more precise, in most of the models, the 
shape parameter, defined by T = fl m h, is typically larger 
than the one in the fiducial model. As a consequence, 
even if we adjust the power spectrum at large scales to 
match the one in the target model, the difference \6Po\ 
can become large as increasing k, leading to a failure of 
the perturbative reconstruction by RegPT-f ast. 

To remedy this situation, a simple but efficient ap- 
proach is to enlarge our set of fiducial models with T- 
parameters that differ from the one of wmap3 model, i.e. 
r = 0.172. Right panel of Fig. [TBI shows the convergence 
results when we supply two extra fiducial models whose 
cosmological parameters are listed in Table ITTl As a fidu- 
cial model with a larger shape parameter, we adopt the 
M001 cosmological model (T = 0.257). Further, for a se- 
cure calculation applicable to general cosmological mod- 
els, we also supply another fiducial model, M023, which 
has a smaller shape parameter (r = 0.139). The initial 
power spectra of those models are plotted in Fig. [T21 de- 
picted as green (M001) and magenta (M023) solid lines. As 
a result, the convergence of the power spectrum calcula- 
tions is dramatically improved, and the RegPT-f ast now 
coincides with rigorous RegPT calculation with < 0.4% 
precision at k < 0.3 /iMpc -1 . Although there still exist 
exceptional cases, in which the fractional difference even- 
tually exceeds 1% precision at k > 0.36 /iMpc -1 , in prac- 
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10 1 1 1 1 — 1 1 — 

io- 2 io-' 

k [h Mpc" 1 ] 

FIG. 12: Linear power spectra Po(k) for 38 cosmological mod- 
els [53l ]. Blue, green, and magenta lines are respectively the 
power spectra of fiducial models wmap3, M001, and M023 used 
for the RegPT-f astcalculation (see Table ILT1 for their cosmo- 
logical parameters). 



tice this is beyond the applicable range of the RegPT cal- 
culation itself. 

With this setting, making use of these three fiducial 
models, RegPT-fast reproduces RegPT direct calcula- 
tions in a wide range of cosmological models and, also 
it does not appear here, for a redshift range of general 
interest. 



B. Comparison with cosmic emulator 

It is now time to discuss the accuracy of the over- 
all RegPT scheme with general cosmic emulator pre- 
dictions. Figs. [14] and [15] summarize the results of the 
comparison for all 38 models, where we plot the ratios 
of power spectra, P(fc)/P no - w i gg i e (fc), at specific redshift 
z = 1. In each panel, magenta solid and black dashed 
lines represent the results of RegPT-fast and the power 
spectrum emulator code, respectively. Also, the fidu- 
cial model used for the RegPT-fast calculation is indi- 
cated, together with the label of the cosmological model. 
The two results mostly coincide with each other, and are 
hardly distinguishable at k < 0.2 /iMpc -1 , where the lin- 
ear theory prediction typically produces a 10% error. At 
k > 0.2 /iMpc -1 , the RegPT-fast results tend to deviate 
from the predictions of the emulator code which proba- 
bly indicates the the limitation of PT treatment. How- 
ever, some models still show a remarkable agreement at 
k < 0.3 h Mpc- 1 (e.g., M009 and M013). 

As the range of applicability of the RegPT scheme de- 
pend on both k and the power spectrum amplitude, fol- 
lowing Refs. (22I [55| . we propose here a phenomenologi- 
cal rule for the domain of applicability of the RegPT cal- 
culations. The proposed upper value for k is fc cr ;t that 
can be obtained from the implicit equation, 

|j jf"* dqP lin (q;z) = C, (42) 

where C is a fixed constant, C = 0.7. The resulting 
values arc depicted as vertical arrows in Figs. [T4l and IT51 
Below the critical wavenumber, the RegPT scheme indeed 
agrees with results of the emulator code, mostly within a 
percent-level precision 12 . We have also checked that this 
is also the case for z = 0.5 with this definition of fc C rit.- 



We however noticed that some models exhibit non-negligible dis- 
crepancy between the results of RegPT-fast and the emulator 
codes, even well below fe cr it. One of such is M015, showing a 
broad-band discrepancy over the plotted range. This is some- 
what surprising in the sense that the RegPT-fast result almost 
converges the linear theory prediction at k < 0. 12 h Mpc - 1 , while 
the result of the emulator code is still away from it. To better un- 
derstand the source of the discrepancy, we have ran TV-body sim- 
ulations for the M015 model - cosmological parameters of M015 
model were set as f2 m = 0.2364, ft b = 0.0384, w = -1.281, 
h = 0.7737, n s = 1.0177, and erg = 0.7692. - with the same 
setup as listed in Table U The resulting power spectrum, esti- 
mated from the ensemble of the 8 independent realizations, is 



18 



0.02 
0.015 
0.01 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
1 fiducial model 



Z=l 



wmap3 EE z = 1 




1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
3 fiducial models 



wmap3 
M001 " 
M023 



0.1 0.2 0.3 
k [h Mpc" 1 ] 



0.1 0.2 0.3 0.4 
k [h Mpc" 1 ] 



FIG. 13: Convergence between the RegPT-f ast treatment and the direct RegPT calculations for 38 cosmological models. The 
fractional difference of the power spectra between RegPT-fast and RegPT direct calculations, PR eg PTfast (k)/Pne S PT (k) — 1, is 
plotted at z = 1. Left and right panels respectively show the results adopting the one and three fiducial models. 



The RegPT scheme is therefore shown to give a fairly 
accurate prediction for the power spectrum in the weakly 
non-linear regime in the sense given above. RegPT direct 
calculations, or (almost) equivalently, RegPT-fast calcu- 
lations with the three fiducial models we prepared, can be 
applied to a wide range of cosmological models. Though 
we did not discuss it here, we expect the same to be also 
true for the correlation function. Finally we note that 
as the relevant scale of weakly non-linear regime grows 
wider for higher redshifts, the applicability and reliability 
of the RegPT scheme is naturally enhanced. On the other 
hand, the emulation schemes to build up interpolated re- 
sults from large sets of TV-body simulations are generally 
efficient in predicting the power spectrum at non-linear 
scales but are more likely to fail at high-z, since the re- 
quirement for the force resolution in iV-body simulation 
becomes more and more severe. In this respect, pertur- 
bative reconstruction schemes such as RegPT - but this 
would also be the case of mptbreeze - are complemen- 
tary to iV-body based predictions. 



VII. CONCLUSION 

It is needless to say that future cosmological observa- 
tions make the development of cosmological tool aiming 
at accurately predicting the large-scale statistical prop- 



supcrposcd in the panel of M015 (green symbols with errorbars) 
and is shown to faithfully trace the RegPT-fast result up to the 
critical wavenumber. It points to a possible flaw in the power 
spectrum emulator, in estimating the smooth power spectrum 
from the ensemble of simulation results, or constructing the in- 
terpolated result of the simulated power spectra. 



erties of the universe highly desirable. In the first part 
of the present paper, based on a renormalized perturba- 
tion theory (PT), we introduced an explicit computation 
scheme applied to the matter power spectrum and cor- 
relation function in weakly non-linear regime that con- 
sistently includes the PT corrections up to the two-loop 
order. The construction of the full expression for the 
power spectrum is based on the T-expansion, i.e. makes 
use of the multi-point propagators which are properly 
regularized so as to recover their expected resummed be- 
havior at high-fc and to match the standard PT result at 
low-fc. We call this regularized PT treatment RegPT. We 
have shown that the RegPT scheme provides an accurate 
prediction for both the power spectrum and the correla- 
tion function, leading to a percent-level agreement with 
./V-body simulations in the weakly non-linear regime. 

In the second half of the paper, we presented a 
method to accelerate the power spectrum calculations. 
The method utilizes prepared data sets for some spe- 
cific fiducial models from which regularized PT calcula- 
tions can be performed for arbitrary cosmological mod- 
els. The main interest of this method is that the eval- 
uation the residual PT corrections between fiducial and 
target cosmological models can be reduced to mere one- 
dimensional integrals. This enables us to dramatically re- 
duce the computational cost, and even with single- node 
calculation by a laptop computer, the power spectrum 
calculation can be done in a few seconds. We call this 
method RegPT-fast, and we have demonstrated that the 
RegPT-fast treatment can perfectly reproduce the direct 
RegPT calculations that involve several multi-dimensional 
integrals. 

We then investigated the range of applicability of the 
RegPT schemes in a broad class of cosmological mod- 
els. For this purpose, we select 38 cosmological mod- 
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FIG. 14: Ratio of power spectra, P(k) /P no -wiggie(fc), at z = 1 for the cosmological models M000-M017. Solid and dotted lines are 
obtained from the RegPT-f ast and cosmic emulator codes, respectively. The fiducial model used for the RegPT-f ast calculation 
is indicated in each panel. The vertical arrows mean the critical wavenumber fc cr it defined by Eq. ((42}, which roughly gives an 
applicable range of RegPT prediction (see text). 



els, and compared the RegPT predictions - eventually in- 
corporating the accelerated computations - with results 
of a power spectrum emulator code, cosmic emulator. 
We show that with the help of three fiducial models the 
RegPT-f ast calculations give reliable predictions for the 
power spectra over this range of cosmological models 13 . 
We furthermore put forward an empirical criterion (|42p 
that gives a good indication of the applicable range of the 
RegPT scheme in k. The RegPT-f ast treatment, together 
with the direct RegPT calculation, has been implemented 



Our analysis is however restricted to flat ui-CDM models. 



in a fortran code that we publicly release as part of this 
paper. 

Although this paper is focused on precision calcula- 
tions of the matter power spectrum, the RegPT frame- 
work as well as the methodology for accelerated calcula- 
tion can naturally be applied to the power spectrum of 
the velocity divergence and the cross-power spectrum of 
velocity and density fields in a similar way. The analy- 
sis of the velocity power spectrum, together with a de- 
tailed comparison with JV-body simulations, will be pre- 
sented elsewhere. Of particular interest is the applica- 
tion of the RegPT schemes to the redshift-space power 
spectrum or correlation function. In this case, not only 
the velocity and density power spectra, but also the 
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FIG. 15: Same as Fig. M but for the models M019-M037. 



multi-point spectra like bispectrum, arising from the non- 
linear mode coupling, seem to play important roles, and 
should be properly modeled. Significance of the effect 
of multi-point spectra has been recently advocated by 
Refs. [HI, [HE-HH , and there appear physical models that 
account for this. Combination of these models with the 
RegPT schemes would be very important, and we will dis- 
cuss it in a near future. 



Acknowledgments 

A.T. is grateful to Toshiya Namikawa for his help- 
ful comments on the code implentation. S.C. thanks 
Christophe Pichon for fruitful comments and discus- 
sion. This work has been benefited from exchange visits 
supported by a bilateral grant from Ministere Affaires 



Etrangeres et Europeennes in France and Japan Soci- 
ety for the Promotion of Science (JSPS). A.T. acknowl- 
edges support from Institutional Program for Young Re- 
searcher Overseas Visit funded by the JSPS. A.T. is also 
supported in part by a Grant-in-Aid for Scientific Re- 
search from the JSPS (No. 24540257). T. N. is supported 
by a Grant-in-Aid for JSPS Fellows (PD: 22-181) and by 
World Premier International Research Center Initiative 
(WPI Initiative), MEXT, Japan. Numerical computa- 
tions for the present work have been carried out in part 
on Cray XT4 at Center for Computational Astrophysics, 
CfCA, of National Astronomical Observatory of Japan, 
and in part under the Interdisciplinary Computational 
Science Program in Center for Computational Sciences, 
University of Tsukuba. F.B and S.C. are also partly sup- 
ported by the French Programme National de Cosmologie 
et Galaxies. 



21 



Appendix A: Perturbative reconstruction of 
RegPT power spectrum at two-loop order 

In this Appendix, we present the set of perturbative 
expressions that are used for the accelerated power spec- 
trum calculation at two-loop order which is implemented 
in the RegPT-f ast code. 

In a similar manner to the one-loop case described in 



Sec. IV Al we can expand the power spectrum expres- 
sion up to two-loop order around the fiducial cosmolog- 
ical model, and obtain the perturbative expression for 
power spectrum in the target cosmological model. Plug- 
ging Eq. (|34|) into the two-loop expression (|23|) and as- 
suming SPq <C fo,fid, the power spectrum is written like 
([551) . and the correction P corr becomes 



Pcorr[fc, ri, a dMrgct ;SP (k)} = 2f« (k; r,) STgUk; rf) P (k) 



d 3 q 
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d 3 pd 3 q ' 



r«(fc;r/)j 5P (k) 

T£l(q, k-q; V )~ P (\k - q\) SP (q) + r£> (g, k - q; r,) 5T^{q, k-q; r,) P (\k - q\) P (q) 

P (\k-p-q\)Po{p)SP (q). (Al) 



r£l(p,q,k-p- q;ri) 



In the above, the perturbations of regularized propaga- 
tors, St{11 and ST^eg, are described as 

«5T« (k; ri) = e 3 " [(1 + a k ) Sf^tk) + e 2 " £rf_L p (fc) 



x e 



(A2) 



Sr[ll(q, k-q; r,) = e 4 " Sr[_> oop (q, k-q) e~ a - , (A3) 

where we define ak = k 2 a\ targot e 2,, /2. The quantities 
<5r|^] oop are defined by 



^ r i-ioo P ( fc ) = 3 



d 3 q 



FW(q,-q,k)6P (q), (A4) 



(2tt) 3 syn 

^-L P W = 30 J f^lF^ m ( qi ,- qi , q2 ,-q 2 ,k) 

Xi 5 0,fid(?l)^o(?2), (A5) 

<JTi-ioo P ( fe i> fe 2) = 6 J -^-^F^l(q,-q,k u k 2 )5P (q). 

(A6) 

The kernels -F S ym are the symmetrized standard PT ker- 
nel for density field. In the above, the angular inte gral s 
are known to be analytically performed (Refs. [23|, 37 1, 
and Bernardeau et al. in prep.), one may write 



^ r i-loop( fc ) 



dqq 2 



f(q;k)SP (q), 



ST 



(i) 

2-loop 



(k) 



2tt 2 
dqi dq 2 q\ q\ 



J(Qi,<l2;k) 



(2tt 2 ) 2 
x Pa^d(qi)SP (92), 



(A7) 



(A8) 



arSLptfci, fc a ) =[^2 K & fc 2, fcs) 5P (g) (A9) 



with the angle-averaged kernels /, J and A" defined by 



(A10) 



= 3/ ^-F^(q,-q,k), 
J(qi,q 2 ;k) = 15 j i_ 

X ^ S ym(9l>-9l, <?2,-<Z2: (AH) 

fci, fca, *s) = 6 / F s ( y t(q, q, fc l5 fc 2 ). (A12) 



Note that fei + &2 = &3- 



The expression for the correction P corr given above 
contains many integrals involving the perturbed linear 
power spectrum, SPq, and some of these require multi- 
dimensional integrals. However, those multi-dimensional 
integration are separately treated, and can be effectively 
reduced to the one-dimensional integrals as follows, 
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«f_L p (fc) 



^.L^(q,k)5Po{q), 



<5r£Lp( fc ) 



d 3 qd 3 q 
(2ir) 3 



= 2 1 |^MW( g ,fc)5Po( g ), 

(g, fc g; »?)] * P ,fid(|fc - <?|) 5P (g) = e 4 " [(1 + a k f J ^ * (2) (<Z, k) SP (q) 



(A13) 
(A14) 



dq q 2 
2tt 2 



Z^(q,k)5P (q)] exp{-2a fe }, 



(A15) 



d 3 q 

(27T) 



3 rP)( 9) fc - q; r,) 6T£l(q, k-q;rj) P 0M (\k - <?|)Po,fid(<7) 



= e 6?? 



j^R {2) (p,k)sp a ( P ) 



-2a k 



(A16) 



d 3 qd 3 q 
(2tt) 6 



r(|(p, 9 ,fe-p-6/ g;? 7)l 2 ^o,fid(|fe-p-Q|)f , o,fid(p)^o(g) = e 6 " e - 2 ^ / S^g, fc) «JP (g). 



(A17) 



In the above, the kernels of the integrals, L, M, A, F, Z, Q, R, and 5, additionally need to be computed, but we 
only have to evaluate them once for each fiducial cosmological model, 



J(p,q; k)P 0M (p), 



X (2) (q, k) = \ J i d N [Fgl(q, k-q)] 2 P M\fk 2 - 2kq ^ + 
^ (2) (<Z, k) = ^J dfi q F^ m (q, k - q) rS oop (q, fc - q) P 0M (\Jk 2 ~ Uq » q + q 2 ), 
Z (2) (q, k ) = \J d »« ff-Lpte, * - ?)] 2 PaM\Jk 2 -2kq» q + q 2 ), 
Q (2) (P,k) = J ^Fgl(q,k-q)K(p:q,\k-q\,k)Po M (\k-q\)P OM (q), 



R {2) (p,k) 



^p- r i-ioo P (9, fe - q)K{p;q, \k - q\,k)P ofid (\k - q\) P ,fid(?), 



S (3) (<?, *) = ^ ^ ^ / 7^ [i^(p, 9, * - P - 9)] 2 P.,fid(|fc - P - <z|) Po,fid(p) 



(A18) 
(A19) 

(A20) 

(A21) 

(A22) 

(A23) 

(A24) 

(A25) 



with the variable \x q defined by fi q = {k- q)/(kq). 

Note that similar to the one-loop case, the correction 
at two-loop order also possesses a one-parameter degree 
of freedom corresponding to a re-scaling the power spec- 
trum amplitude of fiducial model, Po,fid — > cPo,fid- The 
power spectrum difference 5Pq can then be made small 
securing the rapid convergence of this expansion. 

Finally, the data set of kernel functions given above 
are supplemented in the RegPT code, with 301 x 301 log- 
arithmic arrays in (k,q) space. For specific three fidu- 



cial models (i.e., wmap3, M001, and M023), the data have 
been obtained using the method of Gaussian quadrature 
up to three-dimensional integrals and Monte Carlo tech- 
nique for four-dimensional integral. Together with un- 
perturbed part of the PT corrections, these can be used 
as fast calculations of power spectrum at two-loop order. 
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Appendix B: Code description 

In this Appendix, we present a detailed description 
of the fortran code, RegPT, which computes the power 
spectrum and correlation function of density fields valid 
at weakly non-linear regime of gravitational clustering. 



1. Overview 

The code, RegPT, is compiled with the fortran compil- 
ers, if ort or gf ortran. It computes the power spectrum 
in flat wCDM class models based on the RegPT treatment 
when provided with either of transfer function or matter 
power spectrum. It then gives the multiple-redshift out- 
puts for power spectrum, and optionally provides corre- 
lation function data. We have implemented two major 
options for power spectrum calculations: 

• -fast: Applying the reconstruction method de- 
scribed in Sec. IV Al this option quickly computes 
the power spectrum at two-loop level (typically 
a few seconds), using the pre-computed data set 
of PT kernels for fiducial cosmological models. 
We provide the data set for three fiducial models 
(wmap3, M001, and M023, see Table[n|, and the code 
automatically finds an appropriate fiducial model 
to closely match the result of rigorous PT calcula- 
tion with direct-mode. 

• -direct: With this option, the code first applies 
the fast method, and then follows the regularized 



expression for power spectrum (see Eq.[23] with 
regularized propagators *° directly eval- 

uate the multi-dimensional integrals (it typically 
takes a few minutes). The output results are the 
power spectrum of direct calculation and difference 
of the results between fast and direct method. Fur- 
ther, the code gives the data set of PT diagrams 
necessary for power spectrum calculations, from 
which we can construct the power spectrum. We 
provide a supplemental code, read_stf ile . f , with 
which the power spectrum and correlation function 
can be evaluated from the diagram data set in sev- 
eral PT methods, including the standard PT and 
Lagrangian resummation theory (LRT) HE EH as 
well as RegPT treatment (see Appendix 



BTcl) . 



In addition, the code supports the option, 
-direct Hoop, to compute the power spectrum at 
one-loop order. Although this is based on the direct 
calculation with multi-dimensional integration (see 
Eq. [30] with regularized propagators [31] [32]), the one- 
loop expression involves two-dimensional integrals at 
most, and thus the computational cost is less expensive. 
It is potentially useful for the computation of high-z 
correlation function and power spectrum. 



2. Setup 



The RegPT code is available at 



http : //www-utap .phys . s . u-tokyo . ac . jp/~ataruya/regpt_code .html 



A part of RegPT code uses the library for Monte Carlo in- 
tegration, CUBA Before compiling the codes, users 
should download the library package cuba-1 .5, and cor- 
rectly build the file, libcuba.a, compatible with the ar- 
chitecture of user's platform. This can be done in the 
directory /Cuba-1. 5, and just type ./configure and 
make lib. After placing the library file libcuba.a at 
the directory /RegPT/src, users can use the Makefile 
to create the main executable file, RegPT . exe. Note that 
currently available compilers are intel fortran compiler, 
if ort, and GNU fortran compiler, gf ortran. 



putes power spectrum. Users can specify the options in 
the command line, or using the parameter file (suffix of 
file name should be .ini). Sample of parameter file is 
supplied in the code (see directory /RegPT/example). 



3. Running the code 

For running the code with the command-line options, 
Provided with linear power spectrum or transfer func- a simple example is (assuming the code is placed at the 
tion data, the code runs with a set of options, and com- directory, /RegPT) 
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./RegPT.exe -spectrum -infile matterpower_wmap5.dat -nz 2 0.5 1.0 



In the above example, the code first reads the input 
data file, matterpower_wmap5.dat, which is assumed to 
contain linear power spectrum data consisting of two 
columns, i.e., k and Po(k). By default setting, fast 
mode is chosen, and the output result of power spectrum 
is saved to pk_RegPT.dat. With the option -nz 2 0.5 
1 . 0, the output file contains the power spectrum results 
at two redshifts, z = 0.5 and 1.0 (see Appendix IB 4 al for 
output format). Note that by default, the code adopts 
specific values of cosmological parameters. Making use 
of options, users can change the value of cosmological 
parameters appropriately, consistently with input power 
spectrum (or transfer function) data. 

Here we summarize the available options to run the 
code: 

• Verbose level for output message 

-verbose n: This sets the verbose level for output 
information on the progress of numerical compu- 
tation. The available level n is 1 or 2 (default: 
-verbose 1). 

-noverbose: This option suppresses the message 
while running the code. 

• Input data file 

-infile [file]: Input file name of power spec- 
trum or transfer function data is specified (default: 
-infile matterpower.dat). 

-path [ path to input file ] : This specifies the path 
to the input file (default: -path ./). 

-spectrum: With this option, the code assumes 
that the input file is power spectrum data. The 
data consists of two columns, i.e., wavenumber (in 
units of /iMpc^ 1 ) and matter power spectrum (in 
units of h~ 3 Mpc 3 ) (default: -spectrum). The nor- 
malization of power spectrum amplitude can be 
made with the option -sigma8. 

-transfer: With this option, the code assumes 
that the input file is the transfer function data 
created by CAMB. The data should contain 7 
columns, among which the code uses the first and 
seven columns (wavenumber in units of /iMpc -1 
and matter transfer function). The normalization 
of power spectrum amplitude can be made with ei- 
ther of the option -sigma8 or -samp and -spivot. 



• Specification of cosmological parameters 

-sigma8 erg: This option sets the power spectrum 
normalization by as (default: -sigma8 0.817). 
For as < 0, the code will skip the ag normaliza- 
tion. 

-samp ^4 S : This option sets the amplitude of power 
spectrum at pivot scale fc p i vo t (default: -samp 
2.1e-9). This option is used for normalization of 
transfer function data, and is valid when the option 
-transfer is specified. 

-spivot fcpivot^ This option sets the pivot scale 
of CMB normalization in units of Mpc^ 1 (default: 
-spivot 0.05). This option is used for normaliza- 
tion of transfer function data, and is valid when the 
option -transfer is specified. 

-omegam Q m : This option sets the mass density 
parameter (default: -omegam 0.279). This is used 
to estimate the linear growth factor and to compute 
the smooth reference spectrum, P n o-wiggic(fc)- 

-omegab fly,: This option sets the baryon den- 
sity parameter (default: -omegab . 165*omegam). 
This is used to compute the smooth reference spec- 
trum, -Pno-wigglc(fc)- 

-ns n s : This option sets the scalar spectral index. 
This is used to compute the linear power spectrum 
from the transfer function data (option -transfer 
should be specified), and to compute the smooth 
reference spectrum, P no -wiggie(fc)- 

-w w: This option sets the equation of state for 
dark energy (default: -w -1.0). This is used to 
estimate the linear growth factor. 

-h h: This option sets the Hubble parameter (de- 
fault: -h 0.701). This is used to compute the 
power spectrum from the transfer function data, 
and to compute the smooth reference spectrum, 

Pio- wiggle (fc) ■ 

-camb [output parameter file of camb]: With this 
option, the code reads the CAMB output param- 
eter file, and specifies the cosmological parameters 
(O m , Q h , w, h, n s , A s , fcpi vo t). 

• Calculation mode of RegPT 

-fast: This option adopts the fast method of 
power spectrum calculation to give RegPT results. 
This is default setting. 

-direct: This option first applies the fast method, 
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and then follow the direct method for RegPT calcu- 
lation. 

-direct Hoop: With this option, the code adopts 
direct method to compute the power spectrum at 
one-loop order. 

• Setup of fiducial models for fast- and direct-mode 
calculations 

-datapath [path to data directory]: This option 
specifies the path to the data files used for power 
spectrum calculation with fast and direct meth- 
ods (default: -datapath data/). In the directory 
specified with this option, the data set of kernel 
functions given in Appendix [A] and un-perturbed 
part of power spectrum corrections, as well as the 
matter power spectrum should be stored for three 
fiducial cosmological models (wmap3, M001, M023). 

-fiducial [model]: This option sets the specific 
fiducial model among the three, wmap3, M001, and 
M023 (in default setting, the code automatically se- 
lects an appropriate fiducial model). 

• Output data file 

-xicompute: With this option, the code computes 
the correlation function after power spectrum cal- 
culations, and creates the output file. 

-nz n Z\ ■ ■ ■ z n : This option specifies the output 
redshifts for power spectrum calculations. The in- 
teger n specifies the number of redshifts, and subse- 
quent arguments specify the value of each redshift 
(default: -nz 1 1.0). 

-pkfile [file]: This option sets the output 
file name of power spectrum data (default: 
pk_RegPT.dat). 

-xifile [file]: This option sets the output 
file name of correlation function data (default: 
xi_RegPT.dat). 

-stfile [file]: This option sets the output file 
name of PT diagram data (default: st_PT.dat). 

4. Output file format 



pk_RegPT.dat). The columns of this file include 

k, [data for z\], [data for z<j], • • • , [data for z n ] 

The first column is the wavenumber, while the bracket 
[data for Zj\ represents a set of power spectra at given 
redshift Zi and wavenumber k. Number of the data set 
is specified with the option -nz, and each data contains 

Pno- wiggle (k, z^, P\i n {k,Zi), PR Cg pT (k , z t ) , Err(fc) 

Here, the spectrum P no -wiggic is the smooth reference 
spectrum calculated from the no-wiggle formula of lin- 
ear transfer function in Ref. [41(, Py m is the linearly ex- 
trapolated spectrum, and i"R Cg pT(fc, z^ represents the 
power spectrum based on the RegPT calculations with 
fast and/or direct method (depending on the choice of 
options, -fast, -direct or -directlloop). The last 
column, Err, usually sets to zero, but with the option 
-direct, it gives the difference of the power spectra be- 
tween fast and direct methods. 



b. Correlation function data 

With the option -xicompute, the code also provides 
the output file for correlation function data (default file 
name is xi_RegPT.dat). Similar to the power spectrum 
data, the structure of the data is 

r, [data for zi], [data for z%\, • • • , [data for z n ] 

The first column is the separation, while the bracket 
[data for z.;] represents a set of correlation functions given 
at redshift Zi and separation r, containing two columns: 

&m(r,Zi), ^RcgPT^Zi) 

These are simply obtained from the output results of 
power spectrum based on the expression (|2T))) . Note that 
the range of wavenumber for output power spectrum is 
restricted to the wavenumber coverage of input linear 
spectrum (or transfer function). To get a convergent re- 
sult of correlation functions, users may have to supply the 
input data file with a sufficiently wide range of wavenum- 
ber (e.g., 1CT 3 < k < 10/iMpc" 1 ). 



In what follows, wavenumber k and separation r are 
in units of /iMpc -1 and ft, _1 Mpc, respectively All 
the power spectrum data are assumed to be in units of 
h- 3 Mpc 3 . 

a. Power spectrum data 

By default, RegPT code creates the output file 
for the power spectrum data (default file name is 



c. Diagram data 

When users specifies the -direct option, the code ad- 
ditionally provides a set of PT diagram data necessary 
for power spectrum computation, from which we can con- 
struct the power spectrum at one- and two- loop order. 
The output file (default file name is st_PT.dat) includes 
the following columns: 
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k, P 



, (h\ Pi- (k'\ (h\ T { - 1> (k\ p(2)troc-trco/7,\ p(2)tree-lloop/7\ p(2)lloop-Hoop/i.\ p(3) 

\b\K), r\ m \K), L l-loop\ /! 1 2-loop W J -'corr Wi -'corr W) -'corr W) -'cor; 



-d) 



no- wig; 



trcc-trcc 
corr 



(fc) 



Here, the power spectra -P n o-wiggie and Py ln are basically definition |19j). The functions in the remaining four 
the same data as contained in the power spectrum file, columns, P^ ICC ~ tICC , P c ( or r trcc " llo ° p , pio^ 110013-110013 , and 



but these are the extrapolated data at z — (that is, p(3)tree-tree arg j e g ne( j 

Pin corresponds to Pq). The function r^ loop is the 
two-point propagator of the standard PT expansion (see 



Pi 2 r ) r trec - trcc (fc) = 2 J ^L Fs %(q, k q)Fgl(q, k q)P (q)P (\k q\), 
p(2)t r ce-iioo P(fc) = 4 J -f±F^l(q, k - q)T^ loop (q, k q)P (q)Po(\k - q\), 
p c (2)lloop-lloop (fc) = 2 J J0_ T fl op{qt k _ q )T?_\ oop (q, k q)P (q)P a (\k q\), 

p c (3)tree-tree (fc) = g J fP^Lpi^ q . k p q)F^ m (p, q,k-p- q)P Q (p)Po(q)Po(\k p q\ 



Provided the data set above, the power spectrum can 
be constructed with 



J 



P^ p T (k;v)=e^e-^ 



(Bl) 
(B2) 
(B3) 
(B4) 



{l + a fe + e 2 "r^ oop (fc)} P (k) + e^P^ Mee (k) , (B5) 

{i + «* + y + ^ rS„„p(fc) (i + «*) + ^ r^| oop (fc)} 2 Po(fc) 

+ e 2 " {(1 + a k ) 2 P^ ICC - tICC (k) + e 2 "(l + a k )P^ lee - n ° op (k) + e 4 " P c ( 2 r \ lloo P- lloop (fc)} + e 4 " Pi 3 r \ lloo P- lloop (fc) 



(B6) 



r 



for the RegPT calculation at one- and two-loop order, re- 
spectively. Here, a k is given by a k — k 2 a 2 e 2 ' n /2 with era 
being the dispersion of displacement field (see Eq. [16]). 



Note that the diagram data set can be also used to com- 
pute the power spectrum in the standard PT calculations: 



P?5l p (k;r 1 ) = e 2 *P (k)+< 

^Toop(^) V) = ^Woop(^i V) 



2Po(fc)r2 oop (fc) + pi 2 ) 



trco-tree 
corr 



(k) 



(B7) 



J'" i 



PM {ri-L P (fc)} 2 + P^r°- trcc (k) + p^ cc - lloop (k) + 2P (fc) r£[(k) 



(B8) 



With the supplemental code, reacLstf ile . f , users RegPT and standard PT treatments. The code also pro- 
can easily compute the power spectrum in both 
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vides the power spectrum result for LRT [25, 43]. A brief 
instruction on how to run the code and the output format 
of data is described in the header of the code. 

5. Limitation 

Since the RegPT code is the PT-based calculation code 
valid at weakly non-linear scales, the applicability of the 
output results is restricted to a certain range of wavenum- 
ber in power spectrum. We provide an empirical estimate 



of critical wavenumber fc cr it) below which the RegPT re- 
sults are reliable and their accuracy can reach a percent 
level. This is based on Eq. (|42|) with constant value 
C = 0.7(0.3) for two-loop (one-loop) (see Sec. I7IB]| . 
With the option -verbose 2, the code displays the criti- 
cal wavenumbers at output redshifts. Note that the value 
&crit given here is just a crude estimate, and the actual 
domain of applicability may be somewhat wider or nar- 
rower. Users should use the output results with a great 
care. 
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